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The time dependent, isentropic, quasi-one-dimensional equations of gas 
dynamics and other model equations are considered under the constraint of 
characteristic boundary conditions. Analysis of the time evolution shows how 
different initial data may lead to different steady states and how seemingly 
anaraolous behavior of the solution may be resolved. Numerical experimentation 
using time consistent explicit algorithms verifies the conclusions of the 
analysis. The use of implicit schemes with very large time steps leads to 
erroneous results. 


Research was supported in part by the National Aeronautics and Space 
Administration under NASA Contract No. NAS1—17070 while the second and third 
authors were in residence at ICASE, NASA Langley Research Center, Hampton, 
VA 23665. 


fl/85 - H8XWP 


i 



INTRODUCTION 


Consider a steady, isentropic flow in a dual-throat nozzle with equal 
throat areas, and assume that the flow is chocked; then it is well known [1] 
that the flow between the throats can be either completely subsonic or 
supersonic depending on the initial state of the flow and the path taken to 
reach the steady state. If we experiment numerically with the above problem 
using either the isentropic quasi-one-diraensional gas dynamics equation or 
some "simpler" model equation, then some of the results obtained are rather 
peculiar. 

(1) If the initial data corresponds to sufficiently high supersonic flow (or 
sufficiently low subsonic flow), then the steady state flow obtained 
between the two throats is indeed completely supersonic (subsonic). 

(2) If the initial data are completely supersonic (or subsonic), but below a 
certain level (above a certain level) , then the steady state flow 
contains a shock wave connecting the supersonic branch of the solution 
to the subsonic branch. For the model equations considered, the shock 
corresponds to an isentropic jump, and its location depends on the 
initial data. 

(3) Results (1) and (2) above are observed when time accurate schemes are 
used. However, the implicit backwards Euler scheme with large time 
steps yields steady states that are not reachable through a time 
accurate path from any class of nontrivial initial conditions. These 
steady states include not only discontinuous solutions (as observed in 
[2]), but also unstable smooth solutions. 
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(4) The numerical treatment of boundary conditions was very important in 
obtaining the proper results. For example, with central space 
differencing one may have a stable algorithm that does not converge in 
time to a steady state if the sonic conditions are invoked in order to 
supply numerical boundary conditions. 

The purpose of this paper is to present our findings, and to provide, where 
possible, a mathematical explanation of the observed behavior, thereby 
removing the apparent peculiarities. We will show that the nonuniqueness 

aspect of the steady state solution is a by-product of the fact that the 
boundary conditions for the evolution equations are prescribed along 
characteristic curves. This is true for the dual throat problems due to the 
sonic conditions imposed at the throats. The model problems were therefore 
chosen to show this behavior. 

In Section 2, we study the model equation 


G-) -uu-u). 


The relevance of this model equation to the quasi-one-dimensional gas 
dynamical equations is somewhat peripheral. However, it is rich in the number 
of possible steady solutions that it admits, including unstable continuous and 
discontinuous solutions. In this section, we set down the proper way to 
formulate the characteristic boundary conditions for first order quasi-linear 
hyperbolic equations. 
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In Section 3, we consider the model equation 

3u 3 /u 2 \ 

+ 3^ J = sln x cos x - 

This model equation has solutions which qualitatively behave as do those of 
the isentropic dual throat nozzle problem. The simplicity of the model, 
however, affords a detailed study of the possibilities for anomalous 
behavior. This model equation will also show us how to quantify such vague 
terms as sufficiently high (or low) supersonic (subsonic) initial conditions 
that were mentioned in (1) and (2) above. These results are summarized in 
Theorems 1 and 2. 

In Section 4, a model scalar equation is developed which has all of the 
interesting physical aspects of the complete isentropic quasi-one-dimensional 
gas dynamic equations governing the dual throat nozzle problem. To develop 
this equation, our guideline was to retain the differential equation 
exhibiting the characteristic boundary condition, and to model the other 
dependent variable by assuming constant total enthalpy during the time 
evolution. By comparing the theoretical results of the model equation to 
numerical calculations for the complete system of equations, this section 
shows that the proposed single equation is indeed a good model of the complete 
system. Here, by the "goodness" of the model we mean that all of the 
important features of the system are retained. 
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2. FIRST EXAMPLE 

Here we consider the scalar hyperbolic partial differential equation: 


3u 

3t 


+ 


h If-) - "<i - ■>• 


0 < x < 1, t > 0, 


u(x,0) = g(x). 


( 2 . 1 ) 


For reasons mentioned in the introduction, and to be discussed in detail in 
Section 4, we are interested in cases that model physical situations in which 
the boundaries are characteristic. In practice, when (2.1) is solved 
numerically as a characteristic boundary value problem, the boundary 
conditions are imposed dynamically as follows: 


u(0,t) 


u(l ,t) 


0 if u(e Q ,t) > 0 (s Q = Ax) 

.unspecified if u(eg,t) _< 0 
0 if u(e 1 ,t) < 0 (e 1 = 1 - Ax) 

.unspecified if u(e^,t) >^0 


(2.2a) 


(2.2b) 


There are two families of continuous steady states satisfying (2.1) and the 
analytical versions of (2.2): 


u - 0 


(2.3) 


n-x 


u = 1 - e 


(0 < n < l) 


(2.4) 


The stability theory of ordinary differential equations applied to the 
characteristic equation du/dt = u(l - u) easily shows that the steady state 
solution u = 0 is unstable. 

There are also weak solutions connecting various branches (different n's) 
of (2.4). These discontinuous solutions are unstable as will be demonstrated 
now: let 

n l -x 

u L = 1 - e e (2.5) 

be a steady state corresponding to n = , 

n 2 

u R = i - e e x (2.6) 


be another branch. 


Since we want to rule out "expansion shock," i.e., discontinuities that do not 
obey the 'entropy law', we will consider only the case of 1 > r\^ > hj > 0, 
although the analysis is unchanged if ^ . For a steady state shock we 

require u^(xg) + u j^( x g ) = 0. This determines the shock location, xg, to be 

n l . n 2 

„ _ «_ e + e 


We now aslc > what will be the shock speed, x g » y (u L + Ur), if x g is 
perturbed to x s + e? Upon substituting the perturbed shock position in 
(2.5) and (2.6), we get for the new shock speed 


U L + U R 


1 - e~ e « e + 0(e 2 ). 


2 


( 2 . 8 ) 
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Thus, if e > 0 (e < 0) the shock will move to the right (left), showing that 
the solution with a shock is not stable. 

We have thus shown that in the steady state we need consider only the 
smooth solutions in (2.4). We will now demonstrate that these solutions are 
reachable from initial data. The demonstration is first done for the case 
n = 0, g(x) > 0 for all x > 0, and g(0) = 0: 

Consider the problem, (2.1), and let 

g(x) = b( 1 - e x ), b > 0. (2.9) 

The solution to this problem is readily verified as 

u(x,t) = b — — — -^r— . (2.10) 

e _t + b(l - e ) 

Clearly, as t + “, u(x,t) -*■ 1 - e _x , which is a proper steady state. 

Suppose now g(x) is not a multiple of the steady state but is a general 
initial condition still satisfying g(0) = 0, g(x) > 0. The characteristic 
equations are 

4^ = u (2.11) 

dt 

4r = u(l - u) (2.12) 


From (2.12) one gets 


g(Q 

gU) + (1 - g(£) )e -t 


u 


(2.13) 
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where 5 = 5(x,t) is the origin of the characteristic passing through x 
and t. By inserting (2.13) in (2.11) and integrating again along the 
characteristic, we get the following implicit relation between £, x and t: 

e x -5-t = [ g(c) + (i _ g(C))e _t ] (2.14) 


or, upon rearranging 


g(5) = 



(2.15) 


The argument is now as follows: x - 5 is finite (0 x-£ < 1), and thus as 
t g(5) 0, but g(C) 0 only for 5 * 0. Hence, for any finite x, 
as t increases, g(£) takes the large time asymptotic form of 


g(g) ~ ~ 1 (t » 1). (2.16) 

e - 1 

Substituting (2*16) in (2.13) we get 

u(x,t)~- (t » 1). (2.17) 

1 - e C 

Thus, as t + °°, u(x, t) 1 - e x regardless of the detailed form of the 
initial data. 

For other types of initial data (e.g., g(x) = 0 for some x = Xq) the 
proof is the same with n = x Q and the coordinate x transformed to 
x = x - Xq . 

If g(x) has several simple zeros then the interval 0 x 1 is 
subdivided by the zeros. Their relative locations will determine the proper 
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n. In particular, if g(x) is a periodic function, g(xj) = 0, with 
x j , j = 0, 1 • • • ,N, then 


if 


or 


n " *N = 


1 


(i) sgn g'(0) = sgn g'(l) > 0 
(ii) sgn g'(0) = - sgn g'(l) < 0 


(2.18a) 


and 

n = *N-1 
if 

(i) sgn g'(0) = sgn g'(l) < 0 
or 

(ii) sgn g'(0) = - sgn g'(l) > 0, 


(2.18b) 


(where primes denote differentiation with respect to the independent 
variable). In summary, this example demonstrates the richness of possible 
steady state solutions: 

(1) There is an unstable smooth solution, u = 0; 

(2) there are unstable discontinuous solutions; 

(3) there is a one-parameter family of smooth steady states, 

, n-x 
u = 1 - e 

with the value of the parameter depending only on the initial data, a 
direct consequence of the problem having characteristic boundary values. 
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It is interesting to note that if the right-hand side of equation (2.1) 
is taken to be u(u-l), instead of u(l-u), then there is only one possible 
stable steady state solution satisfying the boundary conditions (2.2), namely 


u = 0. 


Note that this was one of the unstable solutions of the previous case. 


2.1 NUMERICAL RESULTS FOR THE FIRST EXAMPLE 

2.1.1 Explicit Form 

The conservative, upwind, first order scheme of Engquist-Osher (E-0) , 
[3] is used to approximate the hyperbolic system of conservation laws 
represented by 

jt + h OH " h(x>u) (2 - 19) 


where h is a source term. Let u^ represent the discrete value of u at 
t n = nAt and x^ = iAx. The explicit, E-0 scheme for equation (2.19) is, 


n+1 




n 
= U, 


1 At 

2 Ax 


2 O “ 6 i+i^ u i+i) + 2 ^ + 6 i-l^ U i-l^ 


+ h(iAx,u")At 


( 2 . 20 ) 


where the switch function 6^ is defined by 
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6 i = 


n 


n. 


u. 




u^ * 0. 


( 2 . 21 ) 


As usual. At satisfies the Courant-Friedrichs-Lewy condition, 


At < 


Ax 


n 

max u. 
l 


( 2 . 22 ) 


and Ax = L/100, where L is the length of the interval of interest. For the 
explicit E-0 scheme convergence was established according to the criterion 


i n+1 n i , . , ~“3 

max u. - u, < 1. x 10 < 

. 1 i l 1 

i 


(2.23) 


The relation given by equation (2.23) is equivalent to requiring the steady 
state operator of (2.20) to be less than 10 . Figure 1 compares the exact 
and computed steady states for equation (2.1) with initial conditions 


g(x) = - sin 2irx. (2.24) 

Note that the steady state satisfies the condition (2.18bi), and that the 
initial and steady state solution is such that no boundary conditions are 


Note that, because of the first order accuracy of the Enquish-Osher scheme, 
Figure 1 shows a slight discrepancy between the analytic and numerical 
solution. The same problem run with Ax = 1/1000 gives results that, on the 
scale of Figure 1, are indistinguishable from the analytic results. This 
comment holds for all other numerical experiments, where, in order to save 
computer time, we used 100 mesh points. 
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imposed at either end of the interval. The same steady state is also obtained 
with 

g(x) = — x ( x— 1 ) ( x - . (2.25) 

Figure 2 compares the exact and computed steady states for initial conditions 

g(x) = sin 2irx. (2.26) 

The steady result is in agreement with the condition (2.18ai). 


2.1.2 Implicit Form 

The slow convergence to steady state characteristic of explicit schemes 
has stimulated research into various acceleration techniques. One of the most 
promising avenues for acceleration consists of recasting the discrete equation 
in implicit form. If we define the increment in time of u by 


*«. * “f 1 - v 


(2.27) 


then the E-0 scheme in implicit form is 


7 Au_ , + 


i+1 i+1 “ u i+l 


(it + 6 i U i " Ax ) Au i “ I (1 + 6 i-l )u i-l Au i-1 


= -l[i 

2 [2 


^ A i+l^ u i+P + A i^ u i^ ~2 (1 + A i-l^ u i-l^ 2 + h(iAx,u^)Ax, 


n.2 1 


n -.2 


(2.28) 
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where 6 is defined as before by equation (2.21). To obtain equation 

(2.28) , terms of order Au^ and higher are neglected. It is easy to see, by 
comparing equations (2.20) and (2.28),. that the right-hand side of equation 

(2.28) is the steady state operator. For the implicit E-0 scheme convergence 
was established by requiring that the steady state operator be less than 10 ^ 
at all mesh points. 

Figure 3 shows the steady state solution obtained using the implicit 
E-0 scheme with 

g(x) = sin 2 ttx (2.29) 

and using infinite Courant number, = 0). The steady state obtained with 
the implicit form of the scheme corresponds to one of the unstable solutions 
of equation (2.1). The stable solution, for g(x) corresponding to equation 

(2.29) , was shown in Figure 2. The peculiar behavior of the implicit 

algorithm at large Courant numbers is further demonstrated in Figure 4 for 

g(x) = - x(x - l)(x ~ -j) (2.30) 

and infinite Courant number. For this case, the steady state reached by 

(2.28) consists of a combination of stable and unstable steady, piecewise 
solutions of equation (2.1). 
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Figure 1. Exact and computed steady states for equation (2.1) with 
initial conditons (2.24) using a time accurate scheme. 
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O Computed 
Exact, unstable 



Figure 3. Exact and computed unstable steady states for equation (2.1) 
with initial conditions (2.29) using an implicit scheme with 
large Courant number. 
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Figure 4. 


Computed steady state for equation (2.1) with initial 
conditions (2.30) using an implicit scheme with large 
Courant number. 
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3. SECOND EXAMPLE 

We now shift our attention to another advection problem. The steady 
states of this problem will have a completely different nature than of those 
found in the previous example. 

The partial differential equation under consideration is 

3 3 2 

gY + -gj- (,^-J = sin x cos x, 0£x£it, t>0 (3.1) 


u(x,0) = g(x), g(0) = g(ir) = 0. 

with boundary conditions as given by (2.2)* 

Here we have two smooth steady state solutions, 

u + = sin x 

(3.2a) 

u = - sin x. 

There is also an infinite number of possible discontinuous solutions of the 
form 

+ ^ 
u = U X < Xg 

, (3.2b) 

U = U X > Xg 

where Xg, the "shock" location, is an arbitrary point in the interval 
(0,ir). Note that, in the steady state, the "shock" speed u c = (u + + u~)/2 

b 

is zero for any 0 < Xg < tt and, therefore, (3.2b) is a legitimate steady 
state solution. In the above solutions we have already eliminated weak 

solutions that violate the 'entropy condition,' u L > u R . 
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We now ask two questions: 

(i) From what class of initial conditions, if any, can either of the two 
smooth solutions, (3.2a), be reached and 
(ii) Under what circumstances is a steady shock established and can its 
location be predicted? 

Consider first the two questions in the particularly simple case when 

g(x) = 3 sin x, (3.3) 

i.e., the initial data is proportional to a smooth steady state. For 8 > 1, 
Theorem 1 will prove that the steady state is the smooth solution u = u + . 
For 3 < -1, a corollary of Theorem 1 leads to u = u~. 

Theorem 1 : The solution of equation (3.1) with boundary conditions 
(2,2 ) y initial conditions (3.3) and 3 > 1 satisfies 

lim u(x,t) = sin x. 


Proof : The characteristic equations resulting from (3.1) are 


dx 

dt 


u 


du _ du _ dF 
dt U dx dx * 


„ 1.2 

F = y sin x; 


(3.4) 


(3.5) 
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Again using £ = 5(x,t) to designate the origin of a characteristic curve 
passing through (x,t), we integrate (3.5) 

\ u 2 - I g 2 (0 = F(x) - F(?) 
or 

u = ±[2F(x) - 2F(5) + g 2 (c) J 1/2 . (3.6) 


As t * 0, C * x and we have to choose the positive branch of (3.6) because 
0 > 1. Thus, using F = (1/2) sin 2 x, 

u = [sin 2 x + (0 2 - 1) sin 2 ?J 1/2 . (3.7) 


We claim now that for t large enough there is a unique correspondence 
between a point (x,t) and £(x,t). in fact, if a shock wave were to appear 
at a certain time t > 0, it will, because of (3.7), separate two positive 
states. The shock wave will have a positive speed, and consequently will 
propagate out of the domain. Therefore, for t large enough, we may 
substitute (3.7) into (3. A), 


or, 


5 [2F(y) - 2FU) + g Z (Oj 


1 1 / 2 


f dy 

5 [sin 2 y + (0 2 - l)sin 2 f;] 1 ^ 2 


(3.8) 

(3.9) 


For every x < ir, the integrand in (3.9) cannot become singular except at the 
lower limit y = 5, £ +■ 0. Thus, t *■ ® as 5 ♦ 0 and the only possible 
solution for very large time is, from (3.7), 
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u --£ [2F(x) - 2F(5) + g 2 (5)J 1/2 = [2F(x) - 2F(0) + g 2 (0)J 1/2 = 
which completes the proof. 

Corollary: Suppose that 3 in (3.3) satisfies 3 < -1 , then 


sin x, 


lim u(x,t) = - sin x. 


Note that in view of (3.8) the results of Theorem 1 hold for any initial 
conditions g(x), such that g(0) = 0, g(x) > sin x. The corollary is thus 
also extended for any g(x) < - sin x. 

Still continuing with the case of g(x) = 3 sin x, we now consider 

0 < 3 < 1. (3.10) 

Here the steady state will be of the form (3.2b). We will show, however, in 
Theorem 2 that the shock location depends on the initial condition. 


Theorem 2: The solution of equation (3.1) with boundary conditions 

(2.2), initial conditions (3.3) and 0 < 3 < 1 satisfies 


lim u(x,t) 


' u + = sin x; 

< 

l u = -sin x; 


0 < x < x g 
Xg < x < n 


(3.11) 
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where 

x g = ir - sin -1 / 1 - p 2 > 1 . (3.12) 

Proof: From the characteristic equation (3.5), with 0 < 3 < 1, we get 

u(x,t) = ±[sin 2 x - (1 - 3 2 )sin 2 (c(x,t))] 1/2 . (3.13) 

In the interval (x — Xg, Xg), Xg as defined in (3.12), u(x,t) cannot change 
sign because the radical in (3.13) cannot vanish in said interval. Since as 
t 0, u(x,t) is positive, we conclude that 

u(x » t ) = [sin 2 X - (1 - 3 2 )sin 2 (5(x,t))] 1/2 , x-x g < x < x g . (3.14) 

In this interval the first characteristic equation ( 3 • 4 ) becomes 

x 

t " / 7 O =? o T7T (3.15) 

5 [ sin y ~ (1 - B 2 )sin 2 (c(x,t))] 1/ 

since t > 0 we must have 5 < x when x-x < x < x . As t ■*■ «, C(x,t) 
must therefore vanish in the limit. It is thus established that 

lim u(x,t) = sin x, (x-x_ < x < x_). (3.16) 

t*X» & b 

Next consider the interval [0,x-x g ). Formally as t », in this leftmost 
interval, 5(x,t) must converge either to zero or x. However, any 
characteristic passing through (x,t) in the interval [0,x-x o ) cannot 

b 

emanate from any 5 > Xg because this would mean a negative slope, and hence 
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a negative u in the interval (x-x g , Xg); this contradicts (3.16). Having 
established that lim £(x,t) = 0, we notice that formally it is possible for 

t+oo 

a characteristic curve, originating in the interval [0,x-x g ), to start with a 
positive slope (required as t -*• 0) and change slope in the interval. This, 
however, will result in a solution containing a "shock" that violates the 
'entropy condition' u L > u R . We thus have our next intermediate result 


lim u(x,t) = sin x, (0 x < Xg). (3.17) 

t-x» 


It now remains for us to show that in the interval Xg < x x 
must be negative and hence equal to - sin x. 

Consider 

/ udx = y [u 2 (0,t) - u 2 (x,t)J. 
or o z 


the solution 


(3.18) 


From (3.13) it follows that u(0,t) = u(x,t) = 0 (in contrast to the case 
of 8 > 1, see (3.7)). Thus 


x x 

j u(x,t)dx = / u(x,0)dx = 20, (3.19) 

0 0 


or, as t 00 

x c 

S 1 . TT 

/ sin x dx + ^ / u(x,t)dx = 20. 

r\ *- 


Performing the integration we find 


x 

lim j u(x,t)dx = 0 - 1 < 0. 
t+ ~ X S 


(3.20) 
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Consider once again (3.15). We see that as t ♦ °°, 5 must go either to 
zero or it. The limit 5(0,®) 0 will always result in a positive u(x,t) 

which is contradicted by (3.20). Therefore, in (3.13) we must choose the 
negative sign, so 


lim u(x,t) = 

t-H» 


- sin x. 


(x g < x < tt), 


(3.21) 


This completes the proof. 


Corollary: Under the conditions of Theorem 2 except that 


-1 < 0 < 0 


the solution still retains the form of (3.11) except that now 


x g = sin 1 / 1 - 0 2 < j . 


(3.22) 


For arbitrary initial data the general behavior is that described in 
Theorems 1 and 2 and their corollaries, i.e., one can get either solution 
(3.2a) or (3.2b). If a "shock" is present in the steady state, the upper and 
lower bounds for its location are given, for g(x) > 0, as follows: 


it - sin 


1 / sin 2 z - g 2 (z) £ x g < it - sin -1 J\ - (/ g(n)dn) 2 , 


(3.23) 
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9 9 

where z maximizes the expression sin^ x - g z (x). For negative initial data 
the bounds are 


sin 




sin 


- g 2 (z) 


< x s < sin" 1 J\-\ 


(/ g(h)dn)‘‘ 

o 


(3.24) 


The upper bound reflects the "area rule" (see (3.18)). The lower bound is the 
first point where u(x,t) can change sign. For g(x) > 0, the upper bound 
becomes sharp (i.e., equals Xg), if u(ir,t) =0 for all t. 

3.1 NUMERICAL RESULTS FOR THE SECOND EXAMPLE 

3.1.1 Explicit Form 

Equation (3.1) is discretized using the explicit E-0 scheme given by 
equation (2.20). Numerical calculations were performed for initial conditions 
given by 

g(x) = 8 sin x, (3.25) 

where 8 is a free parameter such that 0 < 8 ' <_ 1. The steady state shock 
position as a function of 8 is plotted in Figure 5. The numerical results 
are in excellent agreement with the theoretical prediction given by equation 
(3.12). For any 8 > 1, the steady state obtained was u + given by equation 
(3.2a). 

If one uses an algorithm employing central space differencing (e.g., 
MacCormack's scheme), it is then necessary to supply a numerical boundary 
condition. If the steady state value is used for the boundary condition, then 
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the numerical algorithm, though stable, fails to converge to steady state* 
The reason is clearly due to the fact that the numerical boundary condition 
does not allow for a flux through that boundary. As a consequence we have 
(see (3.19)) 

if 

/ u(x,t)dx = 23 
0 


for all t, while the true steady state, u + , requires 


IT 

lim / u(x,t) = 2. 
t-*» 0 


3.1.2 Implicit Form 


Equation (3.1) is discretized using the implicit E-0 scheme given by 
equation (2.28). Once again, numerical calculations were performed for 
initial conditions given by equation (3.25). Now, an additional free 
parameter is 


100 Ax 
ir At 


(3.26) 


which is a measure of how big At is taken in the numerical calculations. 
The results of these series of calculations are given in Figure 6. As 
indicated in the figure, if "small" At's are taken (e>l/ 2 ), then the 
steady state shock location calculated agrees with the theoretical prediction 
of equation (3.12). However, as At increases, the steady state shock 
position is found to the right of its theoretical location. For sufficiently 
high values of At (small e's), the smooth solution is obtained. 
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4. A MODEL FOR QUASI-ONE-DIMENSIONAL FLUID DYNAMICS 

A characteristic boundary value problem, where boundary conditions are 
of the form (2.2), occurs in a double-throat Laval nozzle 



Figure 7. Sketch of double-throat nozzle 


as shown in Figure 7. It is well known [1] that there are two possible smooth 
steady solutions, with sonic conditions at the throats. Between the throats, 
0 < x < 1, the flow can be either completely subsonic or supersonic, the exact 
Mach number distribution, in each case, being dependent on the nozzle area, 
A(x). 

If one considers the isentropic case only, then the flow may be 
described by the quasi-one-dimensional partial differential equations for the 
Riemann variables. 



2 


Y - 1 


= u + 


C> 


<J> = u - 


C, 
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where u is .the velocity, c = (Yp/p) 1/2 is the speed of sound, and y is 
the ratio of specific heats for ideal gases. The equations are: 

|^+ (u + c) ||+ ucF'(x) = 0, (4.1) 

ft + (u _ c) ucF '( x ) = 0, (4.2) 


where F (x) dF(x)/dx d(£nA(x))/dx. This is a hyperbolic system whose 
time evolution is difficult to describe analytically. We therefore seek a 
model for this system, so that with a single equation the most salient 
features are retained. We will present numerical evidence that analytical 
predictions resulting from this.- model equation agree very well with results 
found by numerical integration of the original system (4.1), (4.2). 

The model is derived using a single assumption, namely that the total 
enthalpy is constant not only at steady state but also during the transient 
phase. The mathematical expression of this assumption is: 


2 2 

♦ + r 





(4.3) 


where 


a 


Y + 1 
4 


(4.4) 


c 0 ls the stagnation sound speed, and c* is the sonic sound speed, i.e., 
c* is the sound speed at a sonic throat. 

We now face the choice of solving (4.3) for either i|> in terms of $, 
or vice versa. This dilemma is resolved by recognizing that our "physical" 
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problem will impose characteristic boundary conditions on (4.2), and we would 
like our model equation to retain this feature. Therefore, (4.2) is the 
relevant equation. Solving for tf>, 


♦ - - 




**4 


T 


(2a - 1)* 


2 j 1/2 


(4.5) 


where the positive branch was chosen in order to satisfy the steady state 
boundary condition at x = 0, i.e., at the first throat, where: 





-*> 




2(1 - a) 

2a - 1 C ** 


(4.6) 


Using (4.5) in (4.2), and defining 


At 

♦ ■ ♦/¥* 


(4.7) 


the equation takes the form 


|i+ A(i) ||= H(+)F'(x) (4.8) 


where 


A(<|>) = <J> + 


JLjul 77^7 


(4.9) 


/ 2a - 1 




i - 2 } 2 - 2 < 1 ..--a> ♦ / 1 - 


* / 1 - ♦ 


(4.10) 


/2a - 1 


x = tc*. 


(4.11) 
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Notice that the time scale, t, is determined by the sonic conditions. 

For the sake of clarity let us first examine the simple case of a = 1 
(Y ~ 3), which corresponds to the flow of products caused by detonating solid 

F(x) = £nA(x) • (4.12) 

A 

4>(0) =0 is 

♦ 2 <x) -l (1 - (4.13) 

and so, as in (3.2a) we have two possible steady states; one is positive 
(supersonic), and the other is negative (subsonic): 

* " 2A<x) J (4 - 14) 

♦' - - 

Bearing in mind the results of the previous sections, we will show that in the 

A ^ A 

time evolution problem, <J> and <J> are reachable from different initial 
conditions. Clearly, (4.14) and (4.15) can be connected by a steady shock — 
and again, because of the symmetry of 4> and 4>~, the steady shock location 
could be anywhere in the interval (0,1). We will show that here too bounds 
on Xg can be found and compare them with results of numerical integration of 
the original system (4.1), (4.2). 


explosives. Equation (4.8) then becomes 


$ + ♦ ||- J d - 2} 2 )F'(x); 


A smooth steady state solution of (4.12) with 
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We will concentrate on the positive branch (4.14), showing that if the 
initial condition is given by 

i(x,0) = g(x) = H + = ^[ -^(x) 1 ] 172 (4.16) 

with 

1 < g 2 < , (4.17) 

max 

A A | 

where is the maximum area in the nozzle, then lim <(>(x,t) = <f> (x). A 

t-*°° 

solution of the second characteristic equation, 


a a 



(4.18) 


is given by 

|l - 2* 2 | = |l - 2g 2 (5(x,T))|A(E(x,T))/A(x), (4.19) 


where as before £(x,x) is the origin of a characteristic curve passing 
through (x,t). Since we have chosen (see (4.17)) g 2 (x) to be smaller 

than 1/2, then it follows from (4.19) that 


<|>(x,S) = ± 


A(x) - A(g)fl ~ 2g 2 (g(x,T)) 
2A(x) 


ll/2 


(4.20) 


where £(x,x) is to be determined from the first characterisitc equation 


T 


= ! x ^-. 
K 4>(y,£) 


(4.21) 
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From (4.16) we see that a positive (negative) 0 will initially select a 
positive (negative) branch of (4.20). By an argument similar to that used in 

A 

Theorem 1, it remains for us to show that <j> thus initiated will not change 
sign while evolving to steady state. This follows immediately from (4.20) if 
we use for g(x) equation (4.16) with 0 > 1. 

Next we consider the discontinuous steady state solution. The initial 
data are now taken so that |g(x)| < <J> , see equation (4.14). A lower bound 
for Xg is found by inquiring about the zeros of (4.20) - the argument is the 
same as in the previous section. This happens when 


A(Xg) = A( z ) ( 1 - 2g 2 (z) ) 


(4.22) 


where, as before, z maximizes the expression A(x)(l - 2g 2 (x)). To find the 

upper bound we have to devise an "area rule" for equation (4.12). Because of 

U 

the structure of the right-hand side of (4.12), it is no longer / <J>(x,x)dx 

0 

which is conserved. To find the appropriate "area rule," we divide both sides 

*2 

of (4.12) by 1 - 2<j> >0. The resulting equation after integration by x 

over the interval may be written as: 


§7 / * n dx — / |^[in(l - 2<f> 2 )]dx = — t- F(x) 

1 p) 7 * J o c\ ° A / ~ 

1-/2 4 > 


1 


= 0. 


s~2 o ox ' ■ n 

Under the usual area rule assumptions, 4>(0,t) = <^( 1 , t ) = 0, we have 


(4.23) 


/ in 7 , . + ^~~ 2 — I- dx = const. (4.24) 

0 1 - 7~2 4 


Therefore, an upper bound for Xg is found from: 
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/ ln UlU£ dx = / *„ 1 + O. «<»> dx . (4.25) 

Xg 1 - f~l Ip 0 1 - f~l g(x) 

* _j_ 

When g(x) _< 3<f> > (3 < 1) we expect, as in the previous example, the upper 

and lower bounds on Xg to coincide. This was indeed verified in numerical 

experiments with a particular area distribution, A(x). 

Recalling that (4.12) is a scalar model equation representing the 
systems (4.1), (4.2), we find it interesting to note that this 2x2 system also 
possesses an area rule, namely: 


J 

0 


in 


1 + 


^£dx + 


i - n. p + 


^ / (p+p)dx = I [(i|» 2 (l ,t) + <j. 2 (l,t)) - (i|> 2 (0,t) + 4> 2 (0 , t ) 3 j • (4.26) 


Under the assumption that ({‘(OjO = $(l,t) = 0; ij>(0,t) = ip(l,t), we have 

/ (ij» + <}))dx = 0. (4.27) 


We can now use this to test the "goodness" of our model by comparing the shock 
location predicted from (4.25) with that of the system, whose solution is 
found numerically. This comparison is carried out in the next section. 

Having concluded the analysis of the a = 1 case, let us now return to 
the more general formulation (4.8). In particular, let us consider the case 
of y = 1.4 (a = .6), corresponding to air. We next show how (4.8) may be 
cast in a form similar to the "decoupled" one in (4.12). Multiply both sides 
of (4.8) by r'(<j>) (r' = dr/d<|>): 


3r 

3t 


+ 


r 


3r 

3x 


■ H - l - * (r) ) F'(x) = K(r)(r 

♦'( r) 


r) (r - r_)F'(x), 


(4.28) 


where 


Mr) = /— l - ( /‘ - l - - 1 rKl ~ ^ Kr * 3. (4 29) 

a-r 2 )(r-|/l-| r 2 )(| r + 5 /7T|7) 


r + = /2 ’ 


r = - 


_3 

10 


(4.30) 


The quantities r_ and r + are the values of r which, in the steady state, 
correspond to Mach numbers of zero and infinity, respectively. For general 
values of y, K(r), r + and r_ are replaced by K(r,a), r + (a), r_(a). 
K(r,a) will have the same structure as in (4.29). 

It is easy to verify that K(r), given by (4.29), is a positive, slowly 
monotonically decreasing function in the relevant range r_ < r < r + . In fact 
K(r_) *» 2K(r + ) = .309. In the case of y = 3, i.e., equation (4.12), r = <J> 

and we have r + = — r_ = 1// 2 and K(r) = constant. It is thus clear that 

the topological behavior of (4.28) is the same as that of (4.12), and the 

arguments carry over. In particular the non-unique smooth steady states 

depend on the initial data in the same fashion with respect to 8. 


4.1 NUMERICAL RESULTS FOR QUASI ONE-DIMENSIONAL EQUATIONS 

Here, we study numerically equations (4.1) and (4.2) for y = 3, namely: 
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H + 

Tt + 


3 

37 



(f 2 - 4» 2 )f'(x) . 


(4.32) 


The area of the dual-throat nozzle is defined by 


A(x) 


(1 ~ d) 2 + (l - d(2x ~ l) 2 ) 2 
2(1 - d)(l - d(2x - l) 2 ) 2 


0 < x < 1, (4.33) 


where d is a parameter related to the maximum area by 


= (1 - d ) 2 + 1 
max 2(1 - d) 


(4.34) 


For the numerical experiments, we have used d = 1/6 which results in 


\iax 


1.1. The steady state Mach number distribution is 


M(x) = A(x) ± 


/ 




(4.35) 


and the steady state solution to (4.31) and (4.32) as a function of the Mach 
number is 


* - (1 + M)/(l + M 2 ) 1/2 (4.36) 

4> = - n (1 - M)/(l + M 2 ) 1/2 . (4.37) 


With the stagnation pressure and density used as reference values, the value 
of is /~T. 
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4.1.1 Explicit Form 

Equations (4.31) and (4.32) are discretized using the explicit E-0 
scheme given by equation (2.20). Numerical calculations were performed with 
initial conditions corresponding to 

Mx,°) = 0 X [ ^-/ J , (4.38) 

which is equivalent to (4.16), and with 

♦(*.0) - ^ 1/2 , (4.39) 

or 

*<*.°> - /6 (l - 0 2 («2i^i))V2 . (4.40) 

The initial conditions given by (4.39) correspond to the exact, steady 
solution for while those given by (4.40) correspond to conditions for 
consistent with (4.38) and constant total enthalpy, (4.5). The steady state 
reached was the name in either case; therefore, the results reported here are 
for calculations with (4.40) only. 

Figure 8 summarizes the numerical results. The figure compares the 
predicted steady state shock position as given by (4.25) for the model 
equation (4.12) and the computed position for the system (4.31) and (4.32). 
As is evident from the figure, the agreement is very good* 
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4.1.2 Implicit Form 

Equations (4.31) and (4.32) are discretized using the implicit E-0 
scheme given by equation (2.28). Equations (4.38) and (4.40) are again used 
as initial conditions. The numerical results are summarized in Figure 9. As 
shown in the figure, the steady state shock position depends on the Courant 
number as measured by the parameter 

e = 100 • (4.41) 

At 

For values of e > 10 the steady state shock position is the same as that 
predicted by the explicit form. For values of e < 10 (large At), the steady 
state shock position bifurcates at certain values of 8. 
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Figure 9. Predicted steady state shock position given by (4.25) for 

equation (4.12) and computed position for system (4.31) and 
(4.32) with initial conditions (4.38) and (4.40) using an 
implicit scheme. 
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CONCLUSIONS 

In this paper we analyzed several model equations for characteristic 
initial boundary value problems and examined numerically these as well as the 
quasi-one-dimensional isentropic Euler equations of gas dynamics. 

We show that because of the characteristic nature of the boundary 
conditions the resulting steady states, whether smooth or discontinuous, 
depend on the initial data. Different initial conditions may yield different 
steady states. We also gave an example (see Section 2) of solution to the 
steady state equation which cannot evolve from the initial data. Thus from 
the point of view of the time-dependent equation, we find there are no non- 
unique steady states. 

Another conclusion that one may draw is that in order to have complete 
confidence in the results, numerical schemes for characteristic initial 
boundary value problems should be time consistent and employ only suitable 
boundary conditions. Thus we have shown that implicit methods, even for 
finite Courant numbers, may yield solutions which are piecewise combinations 
of non-unique solutions of the steady state equations. In fact, such 
numerically implicit algorithms may converge to solutions which also include 
parts of unstable steady states. 
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